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Developing simplified, but accurate, theoretical approaches to treat heat transport on all length 
and time scales is needed to further enable scientific insight and technology innovation. Using 
a simplified form of the Boltzmann transport equation (BTE), originally developed for electron 
transport, we demonstrate how ballistic phonon effects and finite-velocity propagation are easily 
and naturally captured. We show how this approach compares well to the phonon BTE, and readily 
handles a full phonon dispersion and energy-dependent mean-free-path. This study of transient 
heat transport shows i) how fundamental temperature jumps at the contacts depend simply on the 
ballistic thermal resistance, ii) that phonon transport at early times approach the ballistic limit 
in samples of any length, and in) perceived reductions in heat conduction, when ballistic effects 
are present, originate from reductions in temperature gradient. Importantly, this framework can 
be recast exactly as the Cattaneo and hyperbolic heat equations, and we discuss how the key to 
capturing ballistic heat effects is to use the correct physical boundary conditions. 


I. INTRODUCTION 


Transient thermal transport is a problem of great in¬ 
terest from both fundamental and applied perspectives. 
For example, it is an important factor in self-heating in 
small electronic devices, in phase-change memory and in 
heat-assisted magnetic recording [I|. Time/frequency- 
domain thermal reflectance (TDTR/FDTR), make use 
of rapid thermal transients to measure the thermal prop¬ 
erties of materials ii- It is now well known that bal¬ 
listic phonon transport becomes important in structures 
with small feature sizes as well as in large structures un¬ 
der rapid transient conditions Si- While first princi¬ 
ples simulations , and other physically detailed tech¬ 
niques, have contributed to our understanding of the ba¬ 
sic science, there remains a need for simplified thermal 
transport approaches that capture the essential physics 
and that are computationally tractable. A technique 
that can be derived from the phonon Boltzmann trans¬ 
port equation (BTE) with clearly identified simplifica¬ 
tions and that provides reasonable accuracy and excellent 
computational efficiency all the way from the ballistic to 
diffusive limits is described in this paper. 


Traditionally, heat transport has been described by 
Fourier’s law with the heat equation (HE), but an un¬ 
physical implication of this app roach is that phonons 
can travel at infinite speed [^. The hyperbolic heat 
equation (HHE) resolves this issue by adding a term 
to the heat equation that ensures a finite propagation 
velocity [13, [HI- It is generally understood that these 
approaches are valid only when heat transport is diffu¬ 
sive and the characteristic length scales are much larger 
than the phonon MEP Sub-continuum approaches 
such as the BTE [13, [HI 13 -[22 1, molecular-dynamics [Hl - 
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and Monte-Carlo simulations [Hl - [23| . capture quasi- 
ballistic phonon transport as well as retardation effects 
due to finite-velocity propagation. It is generally thought 
that the treatment of ballistic phonons requires going be¬ 
yond the HE or HHE [l3 |. 

Sophisticated detailed modeling has added much to the 
fundamental understanding of ballistic effects in phonon 
transport, but there is also a need for simplified, com¬ 
putationally efficient, but accurate analysis techniques. 
Several such approaches have been proposed [i3,[ill,[i3- 
[H, [nil 11313 • These often follow one (or both) of the 
following strategies: i) use the phonon BTE simplified to 
obtain solutions for one specific problem, ii) use a two- 
channel model to treat ballistic and diffusive phonons 
separately. Our goal in this paper is to contribute to 
this work by introducing a technique that combines the 
physics of the phonon Boltzmann equation with the com¬ 
putational efficiency of diffusive equations. 

In a recent paper, we showed that the McKelvey- 
Shockley flux method, a simple form of the BTE, pro¬ 
vides highly accurate solutions for steady-state thermal 
transport from the ballistic to diffusive limits [13] • In this 
paper, we extend the work in [131 to the transient case 
and demonstrate that it naturally includes ballistic and 
diffusive phonon transport and finite-velocity heat prop¬ 
agation, and thus is applicable on all length and time 
scales. Good agreement with the full phonon BTE is ob¬ 
tained, while requiring substantially less computational 
effort. Interestingly, we prove that our simple BTE can 
be rewritten exactly as the Cattaneo and hyperbolic heat 
equations. Both ballistic effects and the finite propaga¬ 
tion time are treated when the correct boundary condi¬ 
tions are used. 

The paper is organized as follows: Section HIl describes 
the McKelvey-Shockley flux method and extends the 
work in [331 to the transient case. In Section mi we 
demonstrate the technique by treating transient heat 
conduction through a dielectric film. Within this section, 
extensions to realistic phonon dispersions and energy- 
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dependent mean-free-paths are discussed and demon¬ 
strated, and the technique is benchmarked against nu¬ 
merical solutions to the full phonon BTE. Section [l^ de¬ 
scribes the connection between the McKelvey-Shockley 
flux method and well-known traditional heat transport 
equations, and lastly, in Section lYl we summarize our 
findings. 


II. THEORETICAL APPROACH 
A. Phonon flux and density 



FIG. 1: Thermal conductor of length L. The mean-free-path 
for backscattering A controls the scattering between the for¬ 
ward/backward fluxes. By specifying the injected phonon 
fluxes (solid arrows), the McKelvey-Shockley flux equations 
describe the evolution of the fluxes inside the material. 


Our starting point is the McKelvey-Shockley flux 
method (s^, an approach developed to treat elec¬ 
tron transport and recently reformulated to handle 
phonon/heat transport [s^. In this work, we consider 
ID transport along x, with y and z directions extending 
to infinity. This technique categorizes phonons into two 
components, those that are forward moving (vx > 0) and 
backward moving {vx < 0). The governing equations de¬ 
scribe how the phonon flux, the product of phonon den¬ 
sity and phonon velocity, vary in space and time. The 
McKelvey-Shockley flux equations are [1^ : 
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where F~^/~ {x,t,e) are the forward/backward phonon 
fluxes [units: ^phonons m“^s“^eV“^], u+(e) is the aver¬ 
age x-projected phonon velocity, A(e) is the mean-free- 
path for backscattering, and e is the phonon energy. 
Eqns. m-m are coupled through their scattering terms 
(right-hand terms), which describe how phonons scatter 
in/out of each flux type. Scattering is described by A, 
which is the average distance traveled along x before a 
forward-moving (‘4-’) or backward-moving (‘—’) phonon 
is backscattered [37|. uj!" is responsible for finite-velocity 
and retardation effects. We refer readers to [s^ for the 
definitions of A and The McKelvey-Shockley flux 
method can be derived from the BTE; the case of steady- 
state is treated in [H, , and a more general derivation 

will be presented in future work. 

Since Eqns. 0-0 are first order in x and t, we must 
specify one spatial and temporal boundary condition for 
each flux component. The spatial boundary conditions 
correspond to specifying the injected phonon fluxes at the 
boundaries: F+(0+, t, e) at the left end, and F~{L~ ,t, e) 
at the right end, where L is the length of the thermal 
conductor (depicted in Fig. [T|). The temporal bound¬ 
ary conditions, depending on the problem at hand, will 
generally correspond to specifying the forward/backward 
fluxes at some given time t': F+(x, t', e) and F~{x, t', e). 

The directed fluxes F^ are related to the net phonon 


flux F and the phonon density n through the relations: 


F{x,t,e) = F+{x,t,e) - F (x,t,e), 
E+(x,t, e) + F~{x,t,e) 


n(x, t, e) = 


Ux (e) 


( 3 ) 

( 4 ) 


The total net phonon flux and total phonon density 
are obtained by integrating over energy: 



Eqns. 0-0, with the appropriate boundary con¬ 
ditions, describe the spatial and temporal evolution of 
phonons in a material. Next, we show how to relate 
phonon fluxes to quantities relevant for heat transport. 


B. Heat current and density 

The heat current is obtained by simply multiplying the 
phonon flux times the phonon energy: 

I^{x,t,e) = eF^{x,t,e). (7) 

We can then define the net heat current (Jg) and heat 
density (Q) as: 

lQ{x,t,e) = I^{x,t,e) - lQ{x,t,e), (8) 

iX (x, t, e) -I- JA (x, t, e) 

Q{Ft,e)= \ (9) 

vl[e) 

By integrating over energy, and hence all phonons, 
we obtain the total heat current (/g*) and heat density 




lQ{x,t,e) de, 
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By replacing phonon flux by heat current Iq, all 
the equations presented in Section III Al can be used to 
treat heat transport. From this point on, unless deemed 
important, the explicit dependence of quantities on en¬ 
ergy e will be dropped to simplify the presentation. Keep 
in mind that a final integration over energy is required 
to compute the total phonon flux/density and total heat 
current/density. 


C. Temperature 


s(/f) 


s{k) 


h=h 


AT=T^-Tj^ 


‘R 


t 


f=0 


In addition to heat current and heat density, it is com¬ 
mon to calculate temperature profiles. We note, how¬ 
ever, that temperature is an equilibrium quantity, and in 
cases where the phonon population is far from equilib¬ 
rium its definition may become ambigious. In this work, 
we will consider small applied temperature differences 
at the ends of the thermal conductor (AT = Tl — Tr, 
where Trji are the temperatures of the left/right con¬ 
tacts), which ensures that the system is near equilibrium 
and temperature is well defined. Note, however, that 
the McKelvey-Shockley flux approach itself is not strictly 
limited to small applied temperature differences. 

The forward/backward heat currents can be expanded 
as: 

+ ( 12 ) 

where Iq = Iq is the equilibrium heat current re¬ 
sulting from a constant background reference tempera¬ 
ture Tref, and SIq is a small correction resulting from a 
small applied AT. 

Temperature is related to heat density by 6Q = Cy ST, 
where Cy is the volumetric heat capacity. Using this with 
Eq. di) and Eq. m, we find that temperature can be 
expressed as: 


ST'^{x,t)-\-ST {x,t) 

+ Tref, 

(13) 

2 

2SI^{x,t) 


(14) 

Cyvt ’ 



where ST^ is the correction in temperature relative to 
Tref for each phonon component (forward and backward) 
and Cy is the heat capacity at Tref. 

The calculational procedure is as follows. Given 
some contact temperatures Tr/Tr (which can be time- 
dependent), the injected heat fluxes from each con¬ 
tact are computed (the relation between Tr/Tr and 
Iq{Q^)/Iq{L~) is shown later). Using the injected heat 
fluxes as boundary conditions, we calculate the directed 
heat fluxes for all x and t using Eqns. replac¬ 

ing with Iq (it is convenient to solve directly for 
SIq since Iq gq adds a constant). If a temperature 
profile is desired, T{x,t) can be determined from Eqns. 
(P^-lfTT)). 


FIG. 2: (Color online) Applied temperature at the left and 
right contacts, Tl and Tr, versus time. Top insets: area 
plots showing occupation of the forward-moving (blue) and 
backward-moving (red) phonon states injected at the left and 
right ends, respectively, of the thermal conductor. For t < 0, 
Tl ~ Tr leading to equal numbers of injected phonons. For 
t > 0, Tl = Tr + AT, which results in a larger number of 
forward states injected at the left side of the material. 

As we will illustrate in the next section, the McKelvey- 
Shockley flux method captures ballistic and finite- 
velocity propagation effects, and is thus applicable on 
all length and time scales. Later we will find that the 
McKelvey-Shockley equations are exactly equivalent to 
well-known diffusive equations. 

III. EXAMPLE: TRANSIENT HEAT 
CONDUCTION ACROSS A DIELECTRIC FILM 

Here we apply the McKelvey-Shockley flux method to 
an example of transient thermal conduction through a 
dielectric film. We assume the dielectric film of length 
L (in which the electronic contribution to thermal trans¬ 
port can be neglected) is joined by two ideal thermalizing 
left and right contacts. As discussed in [s^, an ideal con¬ 
tact is in thermal equilibrium described by Bose-Einstein 
statistics and is perfectly absorbing/reflectionless, mean¬ 
ing phonons that reach the contact do not backscatter at 
the interface. Note that ideal contacts are assumed for 
this example, but this is not a limitation of the McKelvey- 
Shockley flux method; any specified injected heat cur¬ 
rents at the ends of the thermal conductor can be used 
as the boundary conditions. 

Before time t = 0~ the him is at equilibrium Tj-ef = 
300 K, after t = 0+ the left contact temperature is raised 
by some small AT and the right contact temperature is 
kept hxed, as depicted in Fig. O In the case of equilib¬ 
rium contacts, we can easily relate the contact tempera¬ 
ture to the injected heat current: 

/+(0+, t,e) = e /BE(e,Tr.(f)), (15) 

lQ{L-,t,e)=e^-^^^\^hR{e,TR{t)), (16) 

where D{e) is the phonon density of states, /be is the 
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Bose-Einstein distribution and the factor of two comes 
from half the states being forward and backward mov¬ 
ing. Note that the product n+(e)D(e)/2 = M{e)/h, 
where M{e) is the phonon distribution of modes and h 
is Planck’s constant, which can be efficiently calculated 
from the phonon dispersion [dll - ld^ . 

When Tl = Tr, the injected fluxes are equal and no 
net heat current flows. However, if Tl > Tr, the left 
contact will inject a larger heat current than the right 
contact which will drive a net heat current, as illustrated 
in Fig. [2] Using Eqns. (fT^ - (IT6l) as boundary conditions, 
we can solve the McKelvey-Shockley flux equations by 
replacing with Iq in Eqns. o-@. For a small AT = 
Tj^—Tn, we can expand Eq. (HH) to first order around Tr, 
and solve the McKelvey-Shockley flux equations directly 
for SIq, as defined in Eq. (IT^ . 


A. Diamond film: constant velocity and 
mean-free-path 

In this example, we consider diamond as our dielectric 
material, so we may compare to the results of the full 
phonon BTE for the same problem [Toj |. Here, we use 
A = 60 nm Q , which provides a good fit to the steady- 
state and transient results in Ref. [^. Typically A(e) in 
a bulk material can vary by orders of magnitude, however 
in Ref. d, and thus in this work, an effective single 
energy-independent A is used. 

Given diamond’s large Debye temperature 1860 K) 
0, this material can be modeled at 300 K by assum¬ 
ing a linear phonon dispersion with group velocity Vg = 
12 288 m/s. The distribution of modes in this case can 
be written as M(e) = j(Anfi^Vg) [ilj. When consid¬ 
ering a linear dispersion and an energy-independent A, 
the McKelvey-Shockley flux equations need to be solved 
only once instead of at each energy (see Appendix [A|. 
For numerical details see Appendix |B] 

Fig. [3] presents the normalized temperature profile 
(T{x, t) — Tfi/ (Tr — Tr)) versus normalized position xjL 
for diamond hlms of length L = 0.1 ^m (a), 1 ^m (c) 
and 10/rm (e). Lines are solutions of the McKelvey- 
Shockley flux method and symbols are results of the full 
phonon BTE taken from [l0| |. The temperature profiles 
are plotted at different normalized times r = t/tbaii, 
where tbaii = L/vg is the ballistic transit time for a 
phonon traveling at the group velocity Vg. 

An important feature of the temperature profiles is 
the temperature discontinuities at the ends of the ther¬ 
mal conductor. These temperature jumps at the con- 
tacts (dTp) are signatures of ballistic phonon transport 
[lOLIllL l33| , and do not arise in conventional solutions of 
the heat equation or the hyperbolic heat equation [ISliil 
(the HE and HHE would give normalized temperatures 
of 1 and 0 at the left and right contacts, respectively). 
5Tc decreases as L increases and transport shifts from 
ballistic (A >> L) to diffusive (A << L). We previously 
showed that in steady-state dTc = TAT/2 [^, where 
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FIG. 3: (Color online) Normalized temperature profile 
{T{x,t) — Tr)/{Tl — Tr) and normalized heat current 
Iq {x,t)/lQ^^^ in diamond films versus normalized position 
x/L, where is the steady-state ballistic heat current. 
Three film lengths are considered: L = 0.1 pm (a)-(b), l/rm 
(c)-(d), 10 fim (e)-(f). The results are plotted at different nor¬ 
malized times r = t/tbaii, where tbaii = L/vg is the ballistic 
transit time and Vg is the group velocity. Lines are solutions of 
the McKelvey-Shockley flux method, and symbols correspond 
to the phonon BTE p^ . 


T = A/(A -I- T) is the transmission coefficient describing 
the probability of a phonon traveling from one contact 
to the other (ballistic: T —f 1, diffusive: T —f A/T). 
The agreement of dT between our approach and the full 
phonon BTE, over varying L and r, is good; especially 
considering the enormous reduction in complexity of the 
McKelvey-Shockley flux method. We reemphasize that 
our contacts are reflectionless, thus dT is not the result 
of added phonon scattering at the interface - it is the 
signature of ballistic transport. 

Another important characteristic observed in the 
temperature profiles is a finite-velocity propagation of 
phonons. This is most easily seen in the L = 0.1/rm 
film at short times t = 0.1-1, where a wave-front behav¬ 
ior is clearly observed as phonons travel into the thermal 
conductor. A hnite A makes the wave-front decay as it 
moves further into the material; ‘-I-’ phonons scatter and 
become ’ phonons. For an isotropic phonon dispersion 
we have = Vg/2, thus on average phonons have trav- 
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eled only Lj^ at time r = 1, as shown in Fig. [Sja). The 
average time required by the phonon packet to traverse 
the film transitions from (L/n^) for ballistic transport to 
(L^/Z3ph) for diffusive transport, where Dph = 'c+A/2 is 
the phonon diffusion coefficient 

The differences between the McKelvey-Shockley ap¬ 
proach and the full BTE, observed in Fig. |31 are easily 
understood. In the case of L = 0.1/rm, when transport 
is quasi-ballistic, the BTE shows a smoother tempera¬ 
ture profile than McKelvey-Shockley. The temperature 
profile obtained with the BTE is smoother since it con¬ 
siders phonons with Vx varying from 0 to Vg, however 
only a single, angle-averaged is used with McKelvey- 
Shockley leading to a more abrupt temperature distribu¬ 
tion. A finer discretization in angle could be used, but 
that would just make the approach identical to the full 
BTE. The simple discretization into forward and reverse 
streams eliminates two degrees of freedom (a summation 
over A:-points in the 3D Brillouin zone is replaced by an 
integration over energy) while retaining good accuracy 
under transient conditions and producing steady-state 
results that are identical to the BTE [s^ . 

In Fig. m the normalized heat current lQ{x,t)/ 
versus x/L is presented for L = 0.1 p,m (b), 1 p,m (d) 
and 10p,m (f), where = AT*’®'"AT is the steady- 
state ballistic heat current and = CvVx /2 is 

the ballistic thermal conductance [4l|. Interestingly, at 
early times the heat current is much larger than that 
at steady-state. This happens because at very short 
times the ‘-I-’ phonons injected at a; = 0 have not had 
the time to backscatter, leading to enhanced Iq. The 
transient Iq eventually settles to a constant steady-state 
value, which corresponds to the transmission coefficient 
T = Iq! [1^ . In the absence of internal heat gener¬ 
ation, a position-independent steady-state Iq is guaran¬ 
teed with the McKelvey-Shockley flux method; by sub¬ 
tracting Eq. ([2) from Eq. m and using Eqns. (El)-®, 
we obtain the energy balance equation 


dQ _ dlQ 

df dec 


(17) 


Fig. 0] shows the normalized temperature profile and 
normalized heat current versus normalized time for L = 
0.1/rm (a)-(b), 1/rm (c)-(d) and 10/rm (e)-(f). We plot 
the thermal response at the left and right ends of the 
thermal conductor, x = 0 and x = T, respectively. At 
very short times the x = 0 temperature is always close 
to 1/2, then increases and saturates. Similarly the x = 0 
heat current is close to 1 (the ballistic limit), then de¬ 
creases and saturates. Why does this happen? 

Shortly after t = 0~^, the ‘-I-’ states are filled with 
phonons injected from the left contact, while the ’ 
states are empty. From Eq. m, normalized temper¬ 
ature is simply ((5T+ -|- 5T~)/2 (assuming AT = IK for 
simplicity), where at early enough times at the left side 
(5T+ = 1 and 6T~ = 0, which averages to 1/2. This is, in 
fact, the temperature profile of a ballistic thermal con¬ 
ductor (temperature drop occurs only at the contacts) 





r — f ^ ^ball r t / fball 


FIG. 4: (Color online) Normalized temperature profile 
{T{x,t) — Tr)/{Tl — Tr) and normalized heat current 
lQ{x,t)/I q^^^ in diamond films versus normalized time r = 
t/tbaii, where is the steady-state ballistic heat current, 
tbaii = L/vg is the ballistic transit time and Vg is the group 
velocity. Three film lengths are considered: L = 0.1 pm 
(a)-(b), 1pm (c)-(d), 10pm (e)-(f). Lines are solutions of 
the McKelvey-Shockley flux method, and symbols correspond 
to the phonon BTE (with Iq normalized to reproduce the 
McKelvey-Shockley flux value at r = 10“®) [loj |. 


[ 4 ^ . Having only the ‘-I-’ states filled means that the 
heat current is maximal and will at the ballistic limit, as 
shown in Fig. |4l 

Importantly, a normalized temperature approaching 
1/2 accompanied by a heat current near the ballistic limit 
are observed even when L = 10 pm >> A. This shows 
that, at early enough time, ballistic phonon effects are 
important even in “diffusive” samples. Such phenomena 
is relevant for rapid time-resolved thermal characteriza¬ 
tion techniques, including time- and frequency-domain 
thermoreflectance measurements M- 

In Fig. m we find the left side temperature increases 
simultaneously as the heat current decreases. How are 
these two observations related? It is possible to express 
the temperature jump at the contact, physically originat¬ 
ing from nonequilibrium ballistic effects (filled ‘-I-’ states 
and empty ’ states), in terms of a contact resistance 
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(see Appendix ICl) : 


ST,{0,t) 


RTiQjo^t) 

2 


(18) 


where = ^jCyv'^ is the ballistic ther¬ 

mal resistance, and the factor of two comes from divid¬ 
ing the resistance over the left and right contacts. From 
Eq. (1181) . we find that the temperature discontinuity is 
proportional to the heat current at that point, which is 
exactly what we observe when comparing temperature 
and heat current in Fig. 21 

The results presented in Figs. [3] and |4] show that, for 
this simple model problem, a full treatment of the phonon 
BTE agrees surprisingly well with a very simple phonon 
BTE of the form of Eqns. (HJ-©. 


B. Silicon film: full phonon dispersion and 
energy-dependent mean-free-path 

Up to this point, we have considered a diamond film 
using a linear phonon dispersion and energy-independent 
A. For the vast majority of situations, however, accu¬ 
rate quantitative modeling of phonon transport requires 
full phonon dispersions and an energy-dependent A(e). 
Here, we calculate the transient thermal transport prop¬ 
erties of a silicon film, using the phonon band structure 
of bulk Si extracted from first principles (Fig. [S](a)) 
and an energy-dependent A(e) including boundary, defect 
and phonon-phonon scatterings (Fig. ETb)), calibrated 
to experimental data (details found in [33|). The to¬ 
tal (energy-integrated) heat current is obtained from the 
energy-dependent forward/backward heat fluxes Igie) 
using Eq. m- See Appendix [D] for how to calculate 
the total temperature profile from /^(e). 

Fig. [5]presents (c) the normalized temperature and (d) 
the normalized heat current of a L = 30 nm Si film versus 
xjL at different times t. The symbols correspond to re¬ 
sults of the full phonon BTE taken from Ref. , where 
a full phonon dispersion and energy-dependent MFP was 
also used. The overall agreement between the McKelvey- 
Shockley flux method and the full BTE is quite good; 
particularly given the greatly reduced computational ef¬ 
fort of our approach. 

At short times, McKelvey-Shockley flux shows a slower 
propagation speed since we use a single angle-averaged x- 
projected velocity ^^(e) < Vg{e), where Vg{e) is the group 
velocity. Our temperature discontinuities compare very 
well to those of the full BTE, except at the right side at 
short times (as just discussed). Compared to the case of 
diamond, we observe more features in the temperature 
and heat current profiles, because different phonons can 
travel at different velocities and scatter at different rates 
governed by the energy-dependent parameters (e) and 
A(e). We note in passing, that using the average bulk (A) 
leads to significantly different results compared to using 
the energy-dependent A(e). 



xlL xlL 


FIG. 5: (Color online) (a) Bulk Si phonon dispersion calcu¬ 
lated from first principles. Symbols are measured data [d^ . 
(b) MFP distribution versus energy, including boundary, de¬ 
fect and Umklapp scatterings. Dashed line is the calculated 
average bulk A equal to 151 nm. (c) Normalized tempera¬ 
ture profile (T(x, t) — Tr)I{Tl — Tr) and (d) normalized heat 
current lQ{x,t )/in a L = 30nm Si film versus normal¬ 
ized position x/L, where is the steady-state ballistic heat 
current. Lines are solutions of the McKelvey-Shockley flux 
method, and symbols correspond to the phonon BTE [17l] . 


In summary, the McKelvey-Shockley flux method cap¬ 
tures the general trends and essential physics (i.e. tem¬ 
perature jumps and finite-velocity propagation of heat) 
of the phonon BTE, and can be a simple alternative for 
modeling heat transport. 


IV. CATTANEO EQUATION AND THE 
HYPERBOLIC HEAT EQUATION 

We have shown that the McKelvey-Shockley flux 
method can treat heat transport including ballistic ef¬ 
fects, finite-velocity propagation and can provide reason¬ 
able accuracy to the phonon BTE. Here, we demonstrate 
how the McKelvey-Shockley equations can be recast into 
familiar diffusion equations, without approximation. 

By adding the two McKelvey-Shockley flux equations 
(Eqns. ((I])-(I1]), with replaced by Iq), and using the 
definitions for the net heat current (Eq. (|S])) and heat 
density (Eq. ([5])), we find: 

= -.f, (20) 

where tq = \/{2v'^) is the heat relaxation time, Dph = 
uj!"A/2 is the phonon diffusion coefficient and k = 
UyDph = Cvv'^\l2 is the bulk thermal conductivity 
(this expression for k is equal to the traditional relation 
CvVglj'i, where I is the mean-free-path [s^). Eq. (1^ 
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is simply the Cattaneo equation Bini- Combining the 
Cattaneo equation with the energy balance equation (Eq. 
m). we obtain: 

d^T dT ^ d^T 

dt^ ^ dt ~ dx^’ ^ ^ 

where we assume the material parameters are x- 
independent. Eq. m is the well-known hyperbolic heat 
equation, derived from the McKelvey-Shockley flux equa¬ 
tions without making any assumption on L relative to 
A. This indicates that ballistic phonon effects, which are 
present in the McKelvey-Shockley flux approach, are also 
captured by the HHE. Hence, all the results presented in 
this work are in fact simply solutions of the HHE. 

The key to capturing ballistic transport within the 
HHE is to use the correct physical boundary conditions. 
We previously demonstrated that with the proper bound¬ 
ary conditions Fourier’s law and the HE can treat steady- 
state heat transport on all length scales, and accurately 
reproduce results of the phonon BTE . The fact that 
the first order Boltzmann equation can be rewritten ex¬ 
actly as a second order equation is well known in neutron 
transport l48j|. 

Traditionally, one solves the HHE by using the contact 
temperatures, Tl and Tr, as the boundary conditions. 
As we showed earlier, the injected heat fluxes, /+(0+,t) 
and Iq{L~ ,t), are the correct physical boundary condi¬ 
tions. Specifying one flux component at only one end of 
the thermal conductor (Jg at a: = 0+, and Iq a.t x = L~), 
and not both ends, is key to capturing the nonequilibrium 
nature of ballistic transport. 

Although the McKelvey-Shockley flux equations are 
exactly equivalent to the Cattaneo/hyperbolic heat equa¬ 
tion, we find it most convenient to solve the former in¬ 
stead of the HHE. If one desired to use the HHE (in the 
case of ideal contacts), ballistic contact resistances could 
be introduced to effectively capture 5Tc (Eq. (ITOl) and 
Appendix ICl). 

The HHE has previously been derived from the BTE 
[Tsil by assuming local thermodynamic equilibrium at 
each X (we showed earlier that this is the condition of 
diffusive transport), while we found no such approxima¬ 
tion to be necessary. The ballistic-diffusive equations of 
Chen m based on the BTE also yielded an expression 
similar to the HHE, with additional terms due to the bal¬ 
listic phonons. We find that not only is temperature a 
solution of the HHE, but so are net heat current Iq and 
heat density Q which can replace T in Eq. (I^Tl) . The 
HHE derived here suggests that heat travels at a veloc¬ 
ity of \jDph/TQ = as expected from the McKelvey- 
Shockley flux equations. Traditionally, the HHE is solved 
assuming heat propagates (in an isotropic material) at 
a velocity Vg/\f^. However in the case of McKelvey- 
Shockley (which can be derived from the phonon BTE), 
we find an average propagation velocity oi v'^ = Vg 12. 
Lastly, we note that for an isotropic phonon dispersion, 
we have A = (4/3) I = (4/3) Vg r [^, where I is the mean- 
free-path and r is the phonon scattering time. Using this 


relation in rg = A/(2w+), we find rg = (4/3) r, indicat¬ 
ing that the time required for a phonon to scatter from a 
forward-moving state to a backward-moving state is 4/3 
larger than the scattering time. 


V. SUMMARY 

In conclusion, we have shown how to use the 
McKelvey-Shockley flux method to treat transient heat 
transport. By analyzing the transient response of dia¬ 
mond and silicon films, we demonstrated that this ap¬ 
proach (i) captures ballistic phonon effects, such as tem¬ 
perature jumps at the boundaries, (ii) finite-velocity heat 
propagation, and (Hi) can easily support full phonon dis¬ 
persions and energy-dependent MFPs for detailed mod¬ 
eling. Our results show surprisingly good agreement 
with the phonon BTE, while requiring substantially 
less computational effort. Interestingly, we found that 
the McKelvey-Shockley flux equations can be rewritten, 
without approximation, as the Cattaneo and hyperbolic 
heat equations (HHE). The key to capturing ballistic ef¬ 
fects in these diffusive equations is to use the injected 
heat fluxes as the boundary conditions. 

Several physical insights are gained: (i) The temper¬ 
ature jumps at the contacts can be simply captured by 
including a interface resistance equal half the ballistic 
thermal resistance (a material property). The origin of 
the temperature jumps is related to the non-equilibrium 
nature of ballistic transport, and is not the result of en¬ 
hanced scattering at the contacts, (ii) We demonstrate 
that at early enough times, phonon transport can ap¬ 
proach the ballistic limit in samples of any length, which 
is relevant for time-resolved thermal experiments. (Hi) 
We show that the HHE is applicable on all length and 
time scales, which provides an explanation for why dif¬ 
fusion equations work well in situations where ballistic 
effects are present. 

The approach presented in this work is conditioned by 
two approximations. We assume a small applied tem¬ 
perature difference, and we use an angle-averaged x- 
projected phonon velocity at each energy. The latter 
results in wave-front transient heat propagation, rather 
than a smooth one expected in practice due to the wide 
range of Vx from phonons traveling at all angles. 

Given that this approach is strongly connected to fun¬ 
damental physics (it is just a simple BTE derived directly 
from the full BTE), it represents a promising framework 
for analyzing experiments and modeling structures on all 
length and time scales. Finally, this work suggests that 
existing heat transfer simulation tools, based on tradi¬ 
tional diffusive equations, might be modified to capture 
ballistic heat transport. 
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In the case of an energy-independent and A, the 
McKelvey-Shockley flux equations do not need to be eval¬ 
uated at each energy e, and can be simply solved once 
to extract and thus significantly 

simplifying the computational effort. 


Appendix B: Numerical Solution of the 
McKelvey-Shockley Flux Method 


Appendix A: Gray approximation of 
energy-independent and A 


In the case of an energy-independent (i.e. lin¬ 
ear phonon dispersion) and A, the McKelvey-Shockley 
flux equations can be rewritten directly in terms of the 
energy-integrated total heat current (/q*), heat den¬ 
sity (<5*°*) and temperature By integrating the 

McKelvey-Shockley flux equations (Eqns. ([I])-(l2])) over 
energy, we find that Iq can be directly replaced with 


1 

J r+,tot 

-f 

J T+.tot 

T+,tOt 

-'q 

-f 

7-—,tot 

(Al) 

vt 

dt 

dx 

A 

A ’ 

1 

1 7 -,tot 


1 7 -,tot 

7—l-,tOt 

-'q 

-f 

T—,tOt 

(A2) 

Vx 

dt 


dx 

A 

A ■ 


The boundary conditions become 
ir de and Iq [L- ,t) = Iq{L- ,t, e) de. 

By using the definition of total heat density (Eq. ([9|) , and 
Eq. (HD)), we have (5=*=’*°* = Iq’*^°*‘/v^. Hence, in 

Eqns. (TAD-dAD can be replaced with 


1 dQ+’*°* dQ+’*°* 

vt dt dx 

1 dQ-’*°‘ dQ-’*°* 

vt dt dx 


Q+,tot 

A 

Q+,t0t 


A 




(A3) 

(A4) 


The boundary conditions for total heat density are 
Q+'‘°*(0+,t) = x+/^’*°*(0+,t) and = 

,t). Finally, using the following proper¬ 
ties dQ‘°Vdx = d(<5Q*°t)/dx, dQ‘°Vdt = d(<5Q*°t)/dt, 
Q+ - Q- = SQ+ - SQ- and <5r±4ot^ 

Eqns. (IA3I) - (IA4I) can be transformed into: 


1 d(<5T+’‘°*) d(5T+’*°*) 

vt dt dx 

1 d(,5T-4°‘) d((5r-’‘°t) 

Vx dt dx 


A 

jy-l-.tot 

A 


(AS) 

A ’ 
(A6) 


The coupled partial differential equations, that are the 
McKelvey-Shockley Flux equations (Eqns. dD"©), are 
numerically solved using an explicit marching scheme for 
the spatial and temporal derivatives. Given the bound¬ 
ary conditions, it is natural to use a backward differ¬ 
encing in space for the forward-moving flux (Eq. ([p ) 
and a forward differencing in space for the backward- 
moving flux (Eq. ©). A forward differencing in time 
was chosen with both equations, using a time step At < 
[Ax • A/(Ax -I- A)]/x+ to ensure stability. A spatial grid 
resolution of A(x/L) < 0.002 was found to provide good 
accuracy. 


Appendix C: Temperature jumps and ballistic 
thermal resistance 


The left side temperature jump is defined as STc{0, t) = 
Ti — r(0+,t). Noting that = Tfi + (5r+(0+, t), since 
in the contacts the forward/backward temperatures are 
equal, and using Eq. (USD we find 


<5rc(o,t) 


[(5T+(0+,t) -5r-(o+,t)] 
2 


(Cl) 


The net heat current is JQ(x,t) = 5lQ{x,t) — 
which combined with Eq. (HI gives 

Iq ( x , t) = [ST+ (x, t) - ST- (x, t)] , (C2) 

where A*’®'" = Cyv^ j2 is the ballistic thermal conduc¬ 
tance. Combining Eq. (ED and Eq. (ESD, we obtain: 


(5T,(0,f) = (C3) 

where 7?^^" = 1/is the ballistic thermal resistance. 
Thus, the temperature jump can be effectively modeled 
as a contact resistance equal to the ballistic thermal re¬ 
sistance divided by two (for two contacts). An identical 
expression can be derived for the temperature jump at 
the right side {STc{L,t)). 


Appendix D: Energy-averaged temperature profile 


where = /“C'y(e)de and 5T±4ot ^ 

Cy(e)dT=*=(e) de/Cy*. The boundary conditions 
for temperature are i5T+4ot(o+^f) = and 

6T-’^°\L-,t) = TR-T,,f. 


Once the energy-dependent forward/backward heat 
currents Igie) are calculated by solving the McKelvey- 
Shockley flux equations (Eqns. ©-©), the energy- 
dependent temperature profile T{x, t, e) is obtained using 
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Eqns. (Ell-dll. The total (energy-integrated) tempera¬ 
ture is determined using: 




T{x,t,€)Cv{e) de 

rCv(e)de 


(Dl) 


where C'y(e) = e{2M(e)/hv^(e)}ldfBE(e)/dTj and 
Cy * = Cy (e) de is the total bulk heat capacity. 
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